#Libraries and Data

##Libraries

  library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:stats':
## 
##     ar
  library(factoextra)
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
  library(FactoMineR)
  library(corrplot)
## corrplot 0.92 loaded
  library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
  library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
  library(naniar)
  library(rstatix)
## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter
  library(RColorBrewer)
  library(ggeffects)
## 
## Attaching package: 'ggeffects'
## 
## The following object is masked from 'package:cowplot':
## 
##     get_title
  plotopts<-theme_bw() +theme(axis.text=element_text(size=12),axis.title = element_text(size=15),legend.text = element_text(size=15),legend.title = element_text(size=15),strip.text.x = element_text(size=16))

Data

#Pond Data
  pond_wide<-read.csv('Pond Wide Assembled Data.csv',header=T)
    dim(pond_wide)
## [1] 3944   58
#Land Cover Data
  land_cover_wide<-read.csv('Land Cover Wide Assembled Data.csv',header=T,row.names = 1)
#Spawn Climate  
  spawn_climate<-read.csv('Spawn Climate Assembled Data.csv',header=T,row.names = 1)  
#Winter Climate  
  winter_climate<-read.csv('Winter Climate Assembled Data.csv',header=T,row.names = 1)  
#Tadpole Climate  
  tadpole_climate<-read.csv('Tadpole Climate Assembled Data.csv',header=T,row.names = 1)  
  
  
## Year Coverage 
  pond_wide %>% group_by(SITECODE) %>% reframe(min=min(survey_year),max=max(survey_year))
## # A tibble: 10 × 3
##    SITECODE   min   max
##    <chr>    <int> <int>
##  1 T01       1994  2012
##  2 T02       1994  2014
##  3 T03       1994  2011
##  4 T04       1994  2015
##  5 T05       1994  2015
##  6 T06       1994  2010
##  7 T08       1994  1994
##  8 T09       1995  2001
##  9 T10       2001  2015
## 10 T11       1995  2012

Housekeeping

#Variable Labels
  label_df<-data.frame(old=c("year","arable_area","grassland_area","tmax","urban_area","spawn_no3n","forest_area","tmin","af","freshwater_area","other_area","spawn_nh4n","rain"),new=c("Year","Arable Area","Grassland Area","Maximum Temperature","Urban Area","Spawn Nitrate","Forest Area","Minimum Temperature","Air Frost","Freshwater Area","Other Area","Spawn Ammonium","Rain"))

Dataset Assembly

SPAWNDATE data frame (Winter Climate)

#Creating a data frame containing the earliest spawning date, the highest ammonium concentration breeding adults were exposed to???, and the highest nitrate concentration breeding adults were exposed to??? in each ECN location, pond, and year 
  pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% summarize(spawn=mean(spawn_yearday,na.rm=TRUE),
                                                                                                                    prespawn_nh4n=mean(nh4n_numeric[survey_yearday<min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
                                                                                                                 prespawn_no3n=mean(no3n_numeric[survey_yearday<min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
                                                                   ) -> spawndate_df
## Warning: There were 78 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `prespawn_nh4n = mean(...)`.
## ℹ In group 18: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 2011`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 77 remaining warnings.
## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument.
  colnames(spawndate_df)[colnames(spawndate_df) == "survey_year"] <- "year"
  colnames(spawndate_df)[colnames(spawndate_df) == "SITECODE"] <- "location"

#Adding land cover data  
  spawndate_df = merge(x=spawndate_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)
  
#Adding winter climate data   
  spawndate_df = merge(x=spawndate_df,y=winter_climate,by=c("year", "location"),all.x=TRUE)
  
#Replacing infinity values with NA  
  spawndate_df %>% replace_with_na_all(condition = ~.x == -Inf)-> spawndate_df
  spawndate_df %>% replace_with_na_all(condition = ~.x == Inf)-> spawndate_df

HATCHDATE data frame (Spawn Climate)

#Creating a data frame containing the earliest hatching date, the highest ammonium concentration spawn were exposed to???, and the highest nitrate concentration spawn were exposed to??? in each ECN location, pond, and year
  pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% summarize(hatch=mean(hatch_yearday,na.rm=TRUE),
                                                                                                                spawn_nh4n=mean(nh4n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
                                                                                                                 spawn_no3n=mean(no3n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
                                                                  )-> hatchdate_df
## Warning: There were 222 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `spawn_nh4n = mean(...)`.
## ℹ In group 1: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 1994`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.
## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument.
  colnames(hatchdate_df)[colnames(hatchdate_df) == "survey_year"] <- "year"
  colnames(hatchdate_df)[colnames(hatchdate_df) == "SITECODE"] <- "location"
  
#Adding land cover data   
  hatchdate_df = merge(x=hatchdate_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)
  
#Adding data on the climate spawn were exposed to  
  hatchdate_df = merge(x=hatchdate_df,y=spawn_climate,by=c("year", "location"),all.x=TRUE)
  
#Replacing infinity values with NA  
  hatchdate_df %>% replace_with_na_all(condition = ~.x == -Inf)-> hatchdate_df
  hatchdate_df %>% replace_with_na_all(condition = ~.x == Inf)-> hatchdate_df

PERCDEAD data frame

#Creating a data frame containing the highest percentage of spawn found dead, the highest ammonium concentration spawn were exposed to???, and the highest nitrate concentration spawn were exposed to??? in each ECN location, pond, and year
  pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% dplyr::summarize(percdead=mean(percdead_numeric, na.rm=TRUE),
                                                                                                                     spawn_nh4n=mean(nh4n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
                                                                                                                 spawn_no3n=mean(no3n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
                                                                   )-> percdead_df
## Warning: There were 222 warnings in `dplyr::summarize()`.
## The first warning was:
## ℹ In argument: `spawn_nh4n = mean(...)`.
## ℹ In group 1: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 1994`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.
## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument.
  colnames(percdead_df)[colnames(percdead_df) == "survey_year"] <- "year"
  colnames(percdead_df)[colnames(percdead_df) == "SITECODE"] <- "location"

#Adding land cover data
  percdead_df = merge(x=percdead_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)

#Adding data on the climate spawn were exposed to - should I do winter climate instead/as well, a measure of adult fitness?
  percdead_df = merge(x=percdead_df,y=spawn_climate,by=c("year", "location"),all.x=TRUE)

#Replacing infinity values with NA
  percdead_df %>% replace_with_na_all(condition = ~.x == -Inf)-> percdead_df
  percdead_df %>% replace_with_na_all(condition = ~.x == Inf)-> percdead_df

  head(percdead_df)
## # A tibble: 6 × 17
##    year location LCODE percdead spawn_nh4n spawn_no3n forest_area arable_area
##   <int> <chr>    <int>    <dbl>      <dbl>      <dbl>       <dbl>       <dbl>
## 1  1994 T01          1   32.9          NaN        NaN     608706.    9789443.
## 2  1994 T02          1   15            NaN        NaN    1840721.      71881.
## 3  1994 T02          2  NaN            NaN        NaN    1840721.      71881.
## 4  1994 T03          1    0.167        NaN        NaN    1808634.    3341825.
## 5  1994 T04          4    1.43         NaN        NaN      22843.          0 
## 6  1994 T04          1    0            NaN        NaN      22843.          0 
## # ℹ 9 more variables: grassland_area <dbl>, freshwater_area <dbl>,
## #   urban_area <dbl>, other_area <dbl>, tmax <dbl>, tmin <dbl>, af <dbl>,
## #   rain <dbl>, sun <dbl>
  nrow(percdead_df)
## [1] 224
  # 
 ggplot(percdead_df, aes(y=percdead, x=year)) + geom_point(shape=21,size=5,aes(fill=factor(LCODE)))
## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(percdead_df, aes(y=percdead, x=spawn_no3n)) + geom_point(shape=21,size=5,aes(fill=factor(LCODE)))
## Warning: Removed 189 rows containing missing values or values outside the scale range
## (`geom_point()`).

  #  percdead_year

Sample Sizes and Filtering

#HATCH
  hatchdate_df %>% group_by(location) %>%reframe(hatch=length(hatch[which(!is.na(hatch))]))
## # A tibble: 10 × 2
##    location hatch
##    <chr>    <int>
##  1 T01          7
##  2 T02         23
##  3 T03         13
##  4 T04         45
##  5 T05         21
##  6 T06         10
##  7 T08          1
##  8 T09          4
##  9 T10          2
## 10 T11         26
  hatchdate_df<-filter(hatchdate_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & hatch>0)
  nrow(hatchdate_df)
## [1] 145
  table(hatchdate_df$location)
## 
## T01 T02 T03 T04 T05 T06 T11 
##   7  23  13  45  21  10  26
  hatchdate_df$location<-factor(hatchdate_df$location)

#SPAWN  
  spawndate_df %>% group_by(location) %>%reframe(spawn=length(spawn[which(!is.na(spawn))]))
## # A tibble: 10 × 2
##    location spawn
##    <chr>    <int>
##  1 T01         17
##  2 T02         23
##  3 T03         19
##  4 T04         52
##  5 T05         22
##  6 T06         13
##  7 T08          1
##  8 T09          5
##  9 T10          4
## 10 T11         29
  spawndate_df<-filter(spawndate_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & spawn>0)     
  nrow(spawndate_df)
## [1] 175
  spawndate_df$location<-factor(spawndate_df$location)

#PERCDEAD
  percdead_df %>% group_by(location) %>%reframe(percdead=length(percdead[which(!is.na(percdead))]))
## # A tibble: 10 × 2
##    location percdead
##    <chr>       <int>
##  1 T01            16
##  2 T02            22
##  3 T03            14
##  4 T04            67
##  5 T05            18
##  6 T06            13
##  7 T08             0
##  8 T09             5
##  9 T10            14
## 10 T11            28
  percdead_df<-filter(percdead_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & percdead>=0) 
  nrow(percdead_df)
## [1] 178
  table(percdead_df$location)
## 
## T01 T02 T03 T04 T05 T06 T11 
##  16  22  14  67  18  13  28
  percdead_df$location<-factor(percdead_df$location)

PCA and MODELS

SPAWNDATE

Mean Trend Over Time

#Sample Size
 spawndate_df %>% filter(!is.na(spawn)) %>% nrow()
## [1] 175
spawnmeans<- spawndate_df %>% group_by(year,location) %>% reframe(spawn=mean(spawn,na.rm=T))

spawn_siteplot1<-ggplot(spawnmeans,aes(x=year,y=spawn))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + facet_wrap(.~location,nrow=4) + guides(fill="none")  + labs(x="Year",y="Spawn Date (Days 1st Jan)") + plotopts  + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") 
spawn_siteplot1
## `geom_smooth()` using formula = 'y ~ x'

# Correlations
  spawndate_df %>% group_by(location) %>% cor_test(spawn,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2     cor statistic      p conf.low conf.high method   p.adj
##   <fct>    <chr> <chr>  <dbl>     <dbl>  <dbl>    <dbl>     <dbl> <chr>    <dbl>
## 1 T01      spawn year  -0.7      -3.76  0.0019   -0.882   -0.324  Pearson 0.0133
## 2 T02      spawn year  -0.34     -1.63  0.118    -0.657    0.0891 Pearson 0.275 
## 3 T03      spawn year  -0.068    -0.283 0.781    -0.507    0.398  Pearson 0.781 
## 4 T04      spawn year  -0.26     -1.92  0.0605   -0.499    0.0116 Pearson 0.212 
## 5 T05      spawn year   0.28      1.28  0.215    -0.165    0.625  Pearson 0.301 
## 6 T06      spawn year   0.16      0.554 0.59     -0.425    0.656  Pearson 0.688 
## 7 T11      spawn year  -0.25     -1.32  0.199    -0.562    0.133  Pearson 0.301
#LMS
  #Models
   spawn_mods<-spawndate_df %>% group_by(location) %>% do(tidy(lm(spawn~year,data=.))) %>% filter(term=="year") #%>% mutate(p=p.value) %>% mutate(p.adj=p.adjust(p,method="BH"))
  spawn_mods$p.adj<-p.adjust(p=spawn_mods$p.value,"BH")
  
#Annotations  
  ann_text <- data.frame(year = 2005,spawn = 100,lab = "cor = -0.7 \np.adj = 0.01",
                       location = "T01")
spawn_siteplot2<-spawn_siteplot1 + geom_text(data = ann_text,label = ann_text$lab,size=6)
  spawn_siteplot2
## `geom_smooth()` using formula = 'y ~ x'

#Models
  spawn_site_model1<-glm(spawn ~ location*year,data=spawndate_df)
  summary(spawn_site_model1)
## 
## Call:
## glm(formula = spawn ~ location * year, data = spawndate_df)
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       3029.3824  1053.7266   2.875  0.00459 **
## locationT02      -1656.2168  1263.0377  -1.311  0.19162   
## locationT03      -2697.8938  1397.4266  -1.931  0.05529 . 
## locationT04      -1966.8514  1148.6924  -1.712  0.08878 . 
## locationT05      -3765.8058  1274.0553  -2.956  0.00359 **
## locationT06      -3561.9439  1563.8236  -2.278  0.02406 * 
## locationT11      -1909.4359  1355.9569  -1.408  0.16100   
## year                -1.4853     0.5263  -2.822  0.00538 **
## locationT02:year     0.8383     0.6308   1.329  0.18575   
## locationT03:year     1.3530     0.6980   1.938  0.05433 . 
## locationT04:year     1.0042     0.5737   1.751  0.08194 . 
## locationT05:year     1.8642     0.6361   2.930  0.00388 **
## locationT06:year     1.7865     0.7814   2.286  0.02354 * 
## locationT11:year     0.9582     0.6769   1.415  0.15886   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 113.0278)
## 
##     Null deviance: 117574  on 174  degrees of freedom
## Residual deviance:  18197  on 161  degrees of freedom
## AIC: 1339.4
## 
## Number of Fisher Scoring iterations: 2
  plot(ggeffects::ggeffect(spawn_site_model1,term= ~year*location)) 

  confint(spawn_site_model1)
## Waiting for profiling to be done...
##                          2.5 %        97.5 %
## (Intercept)       9.641161e+02  5094.6485864
## locationT02      -4.131725e+03   819.2916000
## locationT03      -5.436800e+03    41.0119295
## locationT04      -4.218247e+03   284.5443839
## locationT05      -6.262908e+03 -1268.7034183
## locationT06      -6.626982e+03  -496.9058782
## locationT11      -4.567063e+03   748.1907367
## year             -2.516893e+00    -0.4536957
## locationT02:year -3.980331e-01     2.0745529
## locationT03:year -1.507074e-02     2.7210279
## locationT04:year -1.201660e-01     2.1286619
## locationT05:year  6.173680e-01     3.1109842
## locationT06:year  2.550412e-01     3.3179075
## locationT11:year -3.685586e-01     2.2848687

Mean Trend with AF Days

spawn_af_plot1<-ggplot(spawndate_df,aes(x=af,y=spawn))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + facet_wrap(.~location,nrow=4) + guides(fill="none")  + labs(x="Air Frost Days",y="Spawn Date (Days 1st Jan)") + plotopts  + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") 
spawn_af_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Correlations
  spawndate_df %>% group_by(location) %>% cor_test(spawn,af,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2      cor statistic      p conf.low conf.high method  p.adj
##   <fct>    <chr> <chr>   <dbl>     <dbl>  <dbl>    <dbl>     <dbl> <chr>   <dbl>
## 1 T01      spawn af    -0.0041   -0.0159 0.987   -0.484      0.477 Pearson 0.987
## 2 T02      spawn af     0.27      1.27   0.217   -0.162      0.612 Pearson 0.401
## 3 T03      spawn af     0.23      0.994  0.334   -0.246      0.622 Pearson 0.468
## 4 T04      spawn af     0.25      1.82   0.0754  -0.0260     0.488 Pearson 0.321
## 5 T05      spawn af     0.4       1.79   0.0917  -0.0689     0.722 Pearson 0.321
## 6 T06      spawn af     0.36      1.27   0.229   -0.240      0.759 Pearson 0.401
## 7 T11      spawn af     0.061     0.320  0.751   -0.312      0.419 Pearson 0.876

Mean Trend with Temperature

spawn_tmax_plot1<-ggplot(spawndate_df)+ geom_smooth(method="lm",aes(x=tmax,y=spawn,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=tmax,y=spawn)) + guides(fill="none")  + labs(x="Mean Winter Temperature",y="Spawn Date (Days 1st Jan)") + plotopts  + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") 
spawn_tmax_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Correlations
  spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2    cor statistic      p conf.low conf.high method   p.adj
##   <fct>    <chr> <chr> <dbl>     <dbl>  <dbl>    <dbl>     <dbl> <chr>    <dbl>
## 1 T01      spawn tmax  -0.31    -1.28  0.219    -0.691    0.196  Pearson 0.307 
## 2 T02      spawn tmax  -0.45    -2.33  0.0298   -0.729   -0.0508 Pearson 0.0695
## 3 T03      spawn tmax  -0.19    -0.804 0.433    -0.594    0.288  Pearson 0.449 
## 4 T04      spawn tmax  -0.34    -2.56  0.0135   -0.561   -0.0746 Pearson 0.0695
## 5 T05      spawn tmax  -0.52    -2.52  0.0218   -0.789   -0.0892 Pearson 0.0695
## 6 T06      spawn tmax  -0.37    -1.31  0.217    -0.764    0.230  Pearson 0.307 
## 7 T11      spawn tmax  -0.15    -0.768 0.449    -0.487    0.233  Pearson 0.449

PCA

#Removing sun
  spawndate_nosun <- spawndate_df %>% dplyr::select(-sun)

#Removing all rows containing NA values
  spawndate_nona <- na.omit(spawndate_nosun)

#Removing response variable
  spawndate_noresp <- spawndate_nona %>% dplyr::select(-spawn)
  
#Removing location variables
  spawndate_for_pca <- spawndate_noresp %>% dplyr::select(-location)
  spawndate_for_pca <- spawndate_for_pca %>% dplyr::select(-LCODE)
  
  head(spawndate_for_pca) 
## # A tibble: 6 × 13
##    year prespawn_nh4n prespawn_no3n forest_area arable_area grassland_area
##   <int>         <dbl>         <dbl>       <dbl>       <dbl>          <dbl>
## 1  1994          0           0.0145    1492094.    5194537.      12178451.
## 2  1996          0.25        0.051       21442.          0       19291786.
## 3  1996          0.04        0.044       21442.          0       19291786.
## 4  1996          0.08        0.0625      21442.          0       19291786.
## 5  1997          0.05        0.02        20741.          0       19292486.
## 6  1997          0.25        0.039       20741.          0       19292486.
## # ℹ 7 more variables: freshwater_area <dbl>, urban_area <dbl>,
## #   other_area <dbl>, tmax <dbl>, tmin <dbl>, af <dbl>, rain <dbl>
#Variable Rename
    spawndatepcalabs<-label_df$new[match(colnames(spawndate_for_pca),label_df$old)]
    spawndatepcalabs[c(2,3)]<-c("Pre-Spawn Ammonium","Pre-Spawn Nitrate")
    colnames(spawndate_for_pca)<-spawndatepcalabs
  
#Performing Principal Components Analysis  
  spawndate_res_pca <- PCA(spawndate_for_pca, graph = FALSE)
  
#Contributions
  dimdesc(spawndate_res_pca, axes = 1:2)
## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Forest Area           0.9629334 1.997195e-13
## Arable Area           0.9577990 7.623573e-13
## Freshwater Area       0.9219618 4.132568e-10
## Urban Area            0.9191159 5.943078e-10
## Other Area            0.8069561 3.284910e-06
## Minimum Temperature   0.6803874 3.535380e-04
## Maximum Temperature   0.6407077 9.885594e-04
## Air Frost            -0.5568433 5.782057e-03
## Grassland Area       -0.9579159 7.408580e-13
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Air Frost             0.6152854 0.0017786579
## Minimum Temperature  -0.5681398 0.0046807595
## Maximum Temperature  -0.6429724 0.0009358000
## Rain                 -0.6650735 0.0005351047
#Producing a rotation plot
  spawndate_rotationplot<-fviz_pca_var(spawndate_res_pca, col.var="black",
   repel = TRUE) +
    labs(x = "PC1", y = "PC2", title = "") + guides(color="none")

Spawndate Data

#Creating a data frame containing the co-ordinates of PC1 and PC2
  spawndate_res_pca$ind$coord
##          Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
## 1   8.26981844 -1.8260942  1.39323957  1.52336216 -1.11175053
## 2  -2.00017040  1.9034884  2.43094661  1.10103625 -0.82950343
## 3  -1.53808478  1.2877144  0.86106097 -0.04125161 -0.18623237
## 4  -1.60423205  1.5265109  1.39354157  0.01797176  0.06829241
## 5  -1.16734789  0.4468061  0.42354312  0.06710644 -0.56053220
## 6  -1.59171770  1.1205533  2.08639707  1.04122984 -0.90227914
## 7  -1.00401547  1.0511366  1.54783192 -0.90571840  1.55804788
## 8   0.37250149 -2.2951325  0.36548375 -0.98074680 -0.47706564
## 9   0.37784815 -2.1896365  0.57738642 -1.07193157 -0.21050733
## 10  0.47516834 -1.4973383  1.92594459 -1.88157243  1.81617688
## 11 -0.15460428 -1.4874470  0.04500407 -0.56480004 -0.88077329
## 12 -0.75125616 -0.6912640 -0.82286615  1.40572963  0.60861209
## 13 -0.71522346 -0.6123265 -0.68593468  1.22876500  0.95086684
## 14 -0.75371127 -0.5230099 -0.47749143  1.29758114  0.98478126
## 15 -0.68515119 -0.5968216 -0.67456220  1.11594801  1.12060558
## 16 -0.70208367 -0.1179275 -0.65409902  0.13905901  0.22346595
## 17  0.02926667 -1.5816522 -0.82418771 -0.63877930 -0.60848361
## 18 -0.81623681  0.7679641 -0.71900411 -1.66413994 -0.55761297
## 19 -0.83966938  0.4861917 -1.27880250 -1.38952464 -1.31037187
## 20 -0.20134918 -1.4209694 -1.85930318  0.65326621  0.46835677
## 21 -1.32980601  1.1104388 -1.98153818  0.16075938 -0.47401854
## 22 -1.34102037  1.1239204 -1.94642650  0.18949414 -0.49299523
## 23  7.67107700  4.0148949 -1.12616402 -0.80284424  0.80292047
  spawndate_pca_data <- spawndate_res_pca$ind$coord[,c(1,2)]

#Creating a data frame containing the earliest spawn date per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
  spawndate_glm_data <- cbind(spawndate_pca_data,spawndate_nona)
  spawndate_glm_data <- as.data.frame(spawndate_glm_data)
  count(spawndate_glm_data) #23 data points
##    n
## 1 23
##  Standardise
  spawndate_glm_data$z_year<-as.numeric(scale(spawndate_glm_data$year))
    spawndate_glm_data$z_tmax<-as.numeric(scale(spawndate_glm_data$tmax))

Correlations

###Model   
  spawndate_glm_data %>% cor_test(.,spawn,Dim.1,method="spearman")
## # A tibble: 1 × 6
##   var1  var2    cor statistic     p method  
##   <chr> <chr> <dbl>     <dbl> <dbl> <chr>   
## 1 spawn Dim.1 -0.45     2936. 0.031 Spearman
  spawndate_glm_data %>% cor_test(.,spawn,Dim.2,method="spearman")
## # A tibble: 1 × 6
##   var1  var2    cor statistic     p method  
##   <chr> <chr> <dbl>     <dbl> <dbl> <chr>   
## 1 spawn Dim.2  0.18     1668. 0.422 Spearman

Multivariate Model

#POND ID
  with(spawndate_df,table(location,LCODE))
##         LCODE
## location  1  2  3  4  5  6
##      T01 17  0  0  0  0  0
##      T02  2 21  0  0  0  0
##      T03 18  1  0  0  0  0
##      T04  1  7  2 10 18 14
##      T05 22  0  0  0  0  0
##      T06 13  0  0  0  0  0
##      T11 16 13  0  0  0  0
  spawndate_df$pond_id<-paste(spawndate_df$location,spawndate_df$LCODE,sep="_")
  spawndate_df_subset<-subset(spawndate_df,!is.nan(spawn) & !is.nan(tmax) & !is.na(tmax))
  nrow(spawndate_df_subset) #172
## [1] 172
  table(spawndate_df_subset$location)
## 
## T01 T02 T03 T04 T05 T06 T11 
##  17  23  19  52  19  13  29
#library(brms)

        bf_spawn<-bf(spawn ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
        bf_spawn_temp<-bf(tmax ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
        #bf_temp_min<-bf(tmin ~ year + (1|p|location)) + gaussian()

###TMAX + TMIN       
  fit_spawn1 <- brm(
  bf_spawn + bf_spawn_temp  + set_rescor(TRUE),
  data = spawndate_df, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99)
)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
  summary(fit_spawn1)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: spawn ~ year + ar(p = 1) + (1 | p | pond_id) 
##          tmax ~ year + ar(p = 1) + (1 | p | pond_id) 
##    Data: spawndate_df (Number of observations: 172) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Correlation Structures:
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_spawn[1]     0.47      0.08     0.32     0.61 1.00     4727     2948
## ar_tmax[1]      0.82      0.05     0.73     0.92 1.00     3258     2431
## 
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 15) 
##                                     Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(spawn_Intercept)                    23.75      4.74    16.41    34.92 1.00
## sd(tmax_Intercept)                      1.70      0.35     1.16     2.55 1.00
## cor(spawn_Intercept,tmax_Intercept)    -0.55      0.18    -0.84    -0.12 1.00
##                                     Bulk_ESS Tail_ESS
## sd(spawn_Intercept)                     1280     2060
## sd(tmax_Intercept)                      1055     2012
## cor(spawn_Intercept,tmax_Intercept)     1277     1641
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## spawn_Intercept   611.37    478.39  -370.24  1508.68 1.00     5209     2897
## tmax_Intercept      9.38     84.72  -155.15   172.16 1.00     3165     2143
## spawn_year         -0.27      0.24    -0.72     0.22 1.00     5226     2897
## tmax_year          -0.00      0.04    -0.08     0.08 1.00     3158     2144
## 
## Further Distributional Parameters:
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_spawn     9.50      0.54     8.52    10.59 1.00     4062     2832
## sigma_tmax      0.64      0.04     0.57     0.72 1.00     4888     3051
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(spawn,tmax)    -0.10      0.08    -0.26     0.06 1.00     4535     2875
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
  #Check 
   pp_check(fit_spawn1,resp='spawn')
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

    pp_check(fit_spawn1,resp='tmax')
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

Multivariate Predictions

fit_spawn1_resid<-resid(fit_spawn1)

spawn_resid_plotdata<- data.frame(cbind(fit_spawn1_resid[,,'spawn'][,c(1,3,4)],fit_spawn1_resid[,,'tmax'][,c(1,3,4)]))
colnames(spawn_resid_plotdata)<-c("spawn","spawn2","spawn97","Temp","Temp2","Temp97")
spawn_resid_plotdata$pond_id<-spawndate_df_subset$pond_id
spawn_resid_plotdata$location<-substr(spawn_resid_plotdata$pond_id,1,3)

spawn_plot1<-ggplot(spawn_resid_plotdata,aes(x=Temp,y=spawn))+ theme_bw() + geom_smooth(method="lm")
spawn_plot2<- spawn_plot1 + geom_errorbar(aes(ymin=spawn2,ymax=spawn97,width=0.15),color="gray77") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=2),color="gray77") + geom_point(data=spawn_resid_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Winter Maximum Temperature Residual",y="Spawn Date Residual",fill="Pond ID") + plotopts  + theme(legend.position = "top") + scale_fill_brewer(palette = "Set2") + annotate("text",label="cor = -0.1 [-0.26,0.05]",x=-2,y=55,size=5)
spawn_plot2
## `geom_smooth()` using formula = 'y ~ x'

#Mean Shift in Spawn (Days) per 1 deg rise
  ggeffects::ggpredict(glm(spawn~Temp,data=spawn_resid_plotdata),terms=list(Temp=1))
## # Predicted values of spawn
## 
## Temp | Predicted |      95% CI
## ------------------------------
##    1 |     -1.58 | -4.18, 1.02

Site Level Posterior Plot

spawn_fit1_site<-ranef(fit_spawn1)$pond_id

spawn_site_plotdata<- data.frame(cbind(spawn_fit1_site[,,'spawn_Intercept'][,c(1,3,4)],spawn_fit1_site[,,'tmax_Intercept'][,c(1,3,4)]))
colnames(spawn_site_plotdata)<-c("spawn","spawn2","spawn97","Temp","Temp2","Temp97")
spawn_site_plotdata$location<-substr(rownames(spawn_site_plotdata),1,3)

spawn_site_plot1<-ggplot(spawn_site_plotdata,aes(x=Temp,y=spawn))+ theme_bw() + geom_smooth(method="lm")  
spawn_site_plot2<- spawn_site_plot1 + geom_errorbar(aes(ymin=spawn2,ymax=spawn97,width=0.15),color="gray55") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=5),color="gray55") + geom_point(data=spawn_site_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Temperature Residual",y="Spawn Date Residual") + plotopts + scale_fill_manual(values=brewer.pal(8,"Set2")[c(1:6,8)]) + annotate("text",label="cor = -0.55 [-0.83,-0.13]",x=-1.5,y=-50,size=6) + theme(legend.position = "top") + labs(fill="")
spawn_site_plot2
## `geom_smooth()` using formula = 'y ~ x'

Graphs

Tests

spawndate_long<- spawndate_df %>% dplyr::select(spawn,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")

  spawndate_corr<-spawndate_long %>% group_by(variable) %>% cor_test(traitval,spawn,method="spearman")
  spawndate_corr<-mutate(spawndate_corr,p.adj=p.adjust(p,"BH"))
  spawndate_corr$Variable<-label_df$new[match(spawndate_corr$variable,label_df$old)]
  spawndate_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% arrange(p.adj)
## # A tibble: 10 × 5
##    Variable              cor statistic        p    p.adj
##    <chr>               <dbl>     <dbl>    <dbl>    <dbl>
##  1 Grassland Area       0.74   185852. 4.27e-29 4.27e-28
##  2 Arable Area         -0.65  1165636. 1.96e-20 9.8 e-20
##  3 Urban Area          -0.62  1446777. 6   e-20 2   e-19
##  4 Maximum Temperature -0.61  1365614. 6.24e-19 1.56e-18
##  5 Minimum Temperature -0.59  1350776. 1.06e-17 2.12e-17
##  6 Other Area          -0.52  1361490. 9.60e-14 1.6 e-13
##  7 Air Frost            0.53   402149. 1.3 e-13 1.86e-13
##  8 Forest Area         -0.53  1083373. 4.62e-13 5.78e-13
##  9 Freshwater Area     -0.49  1333948. 3.97e-12 4.41e-12
## 10 Rain                 0.13   774558. 7.97e- 2 7.97e- 2

Graphs

#Arable area
  spawndate_plot1 <- ggplot(spawndate_df, aes(arable_area/1000000, spawn)) +
    xlab(bquote("Arable area" (km^2))) + ylab("Spawn Date \n(Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.65 ***" ,x=7.5,y=100,size=7)
  spawndate_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Grassland
  spawndate_plot2 <- ggplot(spawndate_df, aes(grassland_area/1000000, spawn)) +
    xlab(bquote("Grassland area" (km^2))) +
    ylab("Spawn Date \n(Days 1st Jan)")  +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2")+ plotopts+ theme(legend.position = "top") + annotate("text",label="cor = 0.74 ***" ,x=8.5,y=100,size=7)
  spawndate_plot2 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Max Temp 
  spawndate_plot3 <- ggplot(spawndate_df, aes(tmax, spawn)) +
    xlab("Mean Maximum Winter Temperature (\u00B0C)") +  ylab("Spawn Date \n(Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = -0.61 ***" ,x=5,y=25,size=7)
  spawndate_plot3
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Air Frost 
  spawndate_plot4 <- ggplot(spawndate_df, aes(af, spawn)) +
    xlab("Air Frost Days") +  ylab("Spawn Date \n(Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.53 ***" ,x=17.5,y=25,size=7)
  spawndate_plot4
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Max Temp Versin 2
  spawndate_plot3_v2 <- ggplot(spawndate_df, aes(tmax, spawn,group=location)) +
    xlab("Mean Maximum Winter Temperature (\u00B0C)") +  ylab("Spawn Date \n(Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = -0.61 ***" ,x=5,y=25,size=7)
  spawndate_plot3_v2
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Air Frost Version 2
  spawndate_plot4 <- ggplot(spawndate_df, aes(af, spawn)) +
    xlab("Air Frost Days") +  ylab("Spawn Date \n(Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.53 ***" ,x=17.5,y=25,size=7)
  spawndate_plot4
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

Combined Plot

Spawn All Data

# spawn_multivariateplot<-plot_grid(spawn_siteplot2,spawn_plot2,nrow=1,labels="AUTO",label_size=25)
# ggsave2('Spawn Multivariate Plot.pdf',spawn_multivariateplot,width=35,height=20,units="cm")

Individual Predictors

#Save PCA Plot
  ggsave2('Spawn PCA PLot.pdf',spawndate_rotationplot,width=10,height=10,units="cm")

#Individual Predictors
  spawndate_combinedplot<-plot_grid(spawndate_plot1,spawndate_plot2,spawndate_plot3,spawndate_plot4,labels=c("C","D","E","F"),nrow=2,label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Stich On Site Specific Temp Trends
  spawndate_combinedplot2<-plot_grid(spawn_siteplot2,spawn_site_plot2,ncol=2,labels = c("A","B"),label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#Stich On Bayesian Model   
  spawndate_combinedplot3<-plot_grid(spawndate_combinedplot2,spawndate_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
  ggsave2('Fig_2_Spawn Grid Plot.pdf',spawndate_combinedplot3,width=30,height=45,units="cm")

HATCHDATE

Mean Trend Over Time

# Correlations
  hatchdate_df %>% group_by(location) %>% cor_test(hatch,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2     cor statistic       p conf.low conf.high method  p.adj
##   <fct>    <chr> <chr>  <dbl>     <dbl>   <dbl>    <dbl>     <dbl> <chr>   <dbl>
## 1 T01      hatch year  -0.23     -0.535 0.616    -0.839     0.631  Pears… 0.677 
## 2 T02      hatch year  -0.36     -1.75  0.0944   -0.671     0.0646 Pears… 0.165 
## 3 T03      hatch year   0.57      2.30  0.0421    0.0274    0.853  Pears… 0.0982
## 4 T04      hatch year  -0.4      -2.87  0.00642  -0.621    -0.121  Pears… 0.0449
## 5 T05      hatch year   0.097     0.423 0.677    -0.350     0.507  Pears… 0.677 
## 6 T06      hatch year   0.45      1.41  0.197    -0.256     0.840  Pears… 0.276 
## 7 T11      hatch year  -0.44     -2.41  0.0239   -0.708    -0.0656 Pears… 0.0836
#Models
   hatch_mods<-hatchdate_df %>% group_by(location) %>% do(tidy(lm(hatch~year,data=.))) %>% filter(term=="year")
   hatch_mods$p.adj<-p.adjust(p=hatch_mods$p.value,"BH")
  hatch_mods
## # A tibble: 7 × 7
## # Groups:   location [7]
##   location term  estimate std.error statistic p.value  p.adj
##   <fct>    <chr>    <dbl>     <dbl>     <dbl>   <dbl>  <dbl>
## 1 T01      year    -0.452     0.845    -0.535 0.616   0.677 
## 2 T02      year    -0.576     0.329    -1.75  0.0944  0.165 
## 3 T03      year     1.32      0.574     2.30  0.0421  0.0982
## 4 T04      year    -0.764     0.267    -2.87  0.00642 0.0449
## 5 T05      year     0.147     0.346     0.423 0.677   0.677 
## 6 T06      year     0.752     0.534     1.41  0.197   0.276 
## 7 T11      year    -1.01      0.419    -2.41  0.0239  0.0835
   # 
    
#Sample Size
 hatchdate_df %>% filter(!is.na(hatch)) %>% nrow()
## [1] 145
hatchmeans<- hatchdate_df %>% group_by(year,location) %>% reframe(hatch=mean(hatch,na.rm=T))

hatch_site_plot1<-ggplot(hatchdate_df,aes(x=year,y=hatch))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~location) + labs(x="Year",y="Hatch Date (Days 1st Jan)") + plotopts  + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") + guides(fill="none")
hatch_site_plot1
## `geom_smooth()` using formula = 'y ~ x'

#Annotations  
  hatch_ann_text <- data.frame(year = 2005,hatch = 60,lab = "cor = -0.4 \np.adj = 0.045",
                       location = "T04")
hatch_site_plot2<-hatch_site_plot1 + geom_text(data = hatch_ann_text,label = hatch_ann_text$lab,size=6)
  hatch_site_plot2
## `geom_smooth()` using formula = 'y ~ x'

Mean Trend with Temperature

hatch_tmax_plot1<-ggplot(hatchdate_df)+ geom_smooth(method="lm",aes(x=tmax,y=hatch,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=tmax,y=hatch)) + guides(fill="none")  + labs(x="Mean Spawning Temperature",y="Hatch Date (Days 1st Jan)") + plotopts  + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") 
hatch_tmax_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Correlations
  spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2    cor statistic      p conf.low conf.high method   p.adj
##   <fct>    <chr> <chr> <dbl>     <dbl>  <dbl>    <dbl>     <dbl> <chr>    <dbl>
## 1 T01      spawn tmax  -0.31    -1.28  0.219    -0.691    0.196  Pearson 0.307 
## 2 T02      spawn tmax  -0.45    -2.33  0.0298   -0.729   -0.0508 Pearson 0.0695
## 3 T03      spawn tmax  -0.19    -0.804 0.433    -0.594    0.288  Pearson 0.449 
## 4 T04      spawn tmax  -0.34    -2.56  0.0135   -0.561   -0.0746 Pearson 0.0695
## 5 T05      spawn tmax  -0.52    -2.52  0.0218   -0.789   -0.0892 Pearson 0.0695
## 6 T06      spawn tmax  -0.37    -1.31  0.217    -0.764    0.230  Pearson 0.307 
## 7 T11      spawn tmax  -0.15    -0.768 0.449    -0.487    0.233  Pearson 0.449

PCA

#Removing sun
  hatchdate_nosun <- hatchdate_df %>% dplyr::select(-sun)

#Removing all rows containing NA values
  hatchdate_nona <- na.omit(hatchdate_nosun)

#Removing response variable
  hatchdate_noresp <- hatchdate_nona %>% dplyr::select(-hatch)

#Removing location variables
  hatchdate_for_pca <- hatchdate_noresp %>% dplyr::select(-location)
  hatchdate_for_pca <- hatchdate_for_pca %>% dplyr::select(-LCODE)
  
#Variable Rename
    hatchdatepcalabs<-label_df$new[match(colnames(hatchdate_for_pca),label_df$old)]
    #hatchdatepcalabs[c(2,3)]<-c("Pre-Spawn Ammonium","Pre-Spawn Nitrate")
    colnames(hatchdate_for_pca)<-hatchdatepcalabs


#Performing Principal Components Analysis     
  hatchdate_res_pca <- PCA( hatchdate_for_pca, graph = FALSE)

#Contributions
  dimdesc(hatchdate_res_pca, axes = 1:2)
## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Arable Area           0.9525272 8.620719e-12
## Urban Area            0.9256264 6.858454e-10
## Maximum Temperature   0.8416292 9.158849e-07
## Minimum Temperature   0.7354440 9.621395e-05
## Spawn Nitrate         0.6012114 3.082123e-03
## Air Frost            -0.6891667 3.889722e-04
## Grassland Area       -0.9254699 6.999493e-10
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Other Area            0.8670282 1.781389e-07
## Freshwater Area       0.8210231 2.842221e-06
## Minimum Temperature   0.6394147 1.355243e-03
## Forest Area          -0.5282852 1.149194e-02
## Air Frost            -0.6774766 5.325217e-04
#Producing a rotation plot   
  hatch_rotationplot<-fviz_pca_var(hatchdate_res_pca, col.var = "black",
                     repel = TRUE
  
                     ) +
    labs(x = "PC1", y = "PC2", title = "") + guides(color="none")

Data

#Creating a data frame containing the co-ordinates of PC1 and PC2
  hatchdate_res_pca$ind$coord
##         Dim.1       Dim.2      Dim.3       Dim.4        Dim.5
## 1   2.1624428  0.25964591 -1.2696770  0.38552765 -0.838354052
## 2  -1.4517261 -0.51039504  0.5089993 -2.09556132 -0.670677701
## 3   2.0207646 -0.45675381 -0.7301879 -0.78668106  0.756124000
## 4  -0.2503202  0.68325087 -1.8213885 -2.02786186  0.437144289
## 5  -0.7179078  0.44545375 -0.9770840 -0.97576641 -1.166821965
## 6   3.0765488  0.30719018 -2.5877631 -1.05572833  2.394646993
## 7  -0.4474291  3.50602841  0.6953934  0.37425726  0.030461439
## 8  -1.0329187  0.06097587 -1.6724709 -0.36589037 -1.573565459
## 9  -1.0421476  0.05873357 -1.6125910 -0.31524971 -1.665392671
## 10 -1.4359662 -1.56826957 -0.4754190  0.06326624  0.888455964
## 11  6.1504394 -2.06732928  3.5165970 -1.48393673 -0.699559162
## 12 -0.3945104  2.78829342  1.1483096 -0.06023258  0.959020762
## 13 -0.4389931  2.77865336  1.4257118  0.17763260  0.531418226
## 14 -0.8021923  2.30953440  1.0896670  0.84135802 -0.108996229
## 15 -0.7959737  2.30895969  1.0693540  0.81834754 -0.073916461
## 16 -1.7241263 -2.07002364  0.1277186  1.08087596  0.393895395
## 17 -2.7973736 -2.98795026  0.8342562  0.73967544  1.212208737
## 18 -2.3445280 -1.50136397  1.0276819 -0.85785737  0.187863205
## 19 -2.3948673 -1.51089401  1.3283574 -0.59602348 -0.278317509
## 20 -1.3346020 -1.52371074 -0.5712024  1.89469024 -0.006008684
## 21  3.2451164 -0.33102203 -1.1280507  2.34477954 -0.606378753
## 22  2.7502704 -0.97900708  0.0737884  1.90037875 -0.103250362
  hatchdate_pca_data <- hatchdate_res_pca$ind$coord[,c(1,2)]

#Creating a data frame containing the earliest hatch date per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
  hatchdate_glm_data <- cbind(hatchdate_pca_data,hatchdate_nona)
  colnames(hatchdate_glm_data) [6]<- c("hatchdate")
  hatchdate_glm_data <- as.data.frame(hatchdate_glm_data)
  count(hatchdate_glm_data) #25 data points
##    n
## 1 22
  table(hatchdate_glm_data$location)
## 
## T01 T02 T03 T04 T05 T06 T11 
##   1   4   0   7   5   0   5
  nrow(hatchdate_glm_data)
## [1] 22

Correlations

###Model   
  hatchdate_glm_data %>% cor_test(.,hatchdate,Dim.1,method="spearman")
## # A tibble: 1 × 6
##   var1      var2    cor statistic        p method  
##   <chr>     <chr> <dbl>     <dbl>    <dbl> <chr>   
## 1 hatchdate Dim.1  -0.7     3009. 0.000293 Spearman
  hatchdate_glm_data %>% cor_test(.,hatchdate,Dim.2,method="spearman")
## # A tibble: 1 × 6
##   var1      var2    cor statistic     p method  
##   <chr>     <chr> <dbl>     <dbl> <dbl> <chr>   
## 1 hatchdate Dim.2 0.033     1712. 0.883 Spearman

Multivariate Model

#POND ID
  with(hatchdate_df,table(location,LCODE))
##         LCODE
## location  1  2  4  5  6  7
##      T01  7  0  0  0  0  0
##      T02  2 21  0  0  0  0
##      T03 12  1  0  0  0  0
##      T04  0  5  8 18 13  1
##      T05 21  0  0  0  0  0
##      T06 10  0  0  0  0  0
##      T11 15 11  0  0  0  0
  hatchdate_df$pond_id<-paste(hatchdate_df$location,hatchdate_df$LCODE,sep="_")
  hatchdate_df_subset<-subset(hatchdate_df,!is.nan(hatch) & !is.nan(tmax) & !is.na(tmax))
  nrow(hatchdate_df_subset) #143
## [1] 143
#library(brms)

        bf_hatch<-bf(hatch ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
        bf_hatch_temp<-bf(tmax ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
        #bf_temp_min<-bf(tmin ~ year + (1|p|location)) + gaussian()

###TMAX + TMIN       
  fit_hatch1 <- brm(
  bf_hatch + bf_hatch_temp  + set_rescor(TRUE),
  data = hatchdate_df_subset, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99)
)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
  summary(fit_hatch1)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: hatch ~ year + ar(p = 1) + (1 | p | pond_id) 
##          tmax ~ year + ar(p = 1) + (1 | p | pond_id) 
##    Data: hatchdate_df_subset (Number of observations: 143) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Correlation Structures:
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_hatch[1]     0.34      0.08     0.18     0.50 1.00     3959     2917
## ar_tmax[1]      0.66      0.07     0.52     0.80 1.00     3854     2810
## 
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 14) 
##                                     Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(hatch_Intercept)                    22.26      4.55    15.28    32.95 1.00
## sd(tmax_Intercept)                      1.49      0.33     1.00     2.27 1.00
## cor(hatch_Intercept,tmax_Intercept)    -0.47      0.21    -0.80     0.02 1.00
##                                     Bulk_ESS Tail_ESS
## sd(hatch_Intercept)                     1421     2085
## sd(tmax_Intercept)                      1286     1864
## cor(hatch_Intercept,tmax_Intercept)     1437     1851
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## hatch_Intercept   758.96    437.35  -111.00  1618.20 1.00     4250     2971
## tmax_Intercept    -44.80     50.75  -144.45    49.17 1.00     3606     2464
## hatch_year         -0.33      0.22    -0.76     0.11 1.00     4256     2989
## tmax_year           0.03      0.03    -0.02     0.08 1.00     3609     2464
## 
## Further Distributional Parameters:
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_hatch     9.82      0.62     8.71    11.12 1.00     3957     2926
## sigma_tmax      0.60      0.04     0.53     0.68 1.00     4258     2895
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(hatch,tmax)    -0.37      0.08    -0.52    -0.20 1.00     3653     2898
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
  bayes_R2(fit_hatch1)
##          Estimate  Est.Error      Q2.5     Q97.5
## R2hatch 0.8409414 0.01221426 0.8126497 0.8600146
## R2tmax  0.8497104 0.01190528 0.8226606 0.8682156
  #Check 
  # pp_check(fit_spawn1,resp='highestpercdeadlogit')
  #  pp_check(fit_spawn1,resp='tmax')

Multivariate Predictions

fit_hatch1_resid<-resid(fit_hatch1)

hatch_resid_plotdata<- data.frame(cbind(fit_hatch1_resid[,,'hatch'][,c(1,3,4)],fit_hatch1_resid[,,'tmax'][,c(1,3,4)]))
colnames(hatch_resid_plotdata)<-c("hatch","hatch2","hatch97","Temp","Temp2","Temp97")
hatch_resid_plotdata$pond_id<-hatchdate_df_subset$pond_id
hatch_resid_plotdata$location<-hatchdate_df_subset$location

hatch_plot1<-ggplot(hatch_resid_plotdata,aes(x=Temp,y=hatch))+ theme_bw() + geom_smooth(method="lm")  
hatch_plot2<- hatch_plot1 + geom_errorbar(aes(ymin=hatch2,ymax=hatch97,width=0.15),color="gray77") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=2),color="gray77") + geom_point(data=hatch_resid_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Spawning \nTemperature Residual",y="Hatch Date Residual",fill="Location") + plotopts + theme(legend.position = "top") + scale_fill_brewer(palette="Set2") + annotate("text",x=-1.8,y=-25,label="cor = -0.37 \n[-0.52, -0.22]",size=6)
hatch_plot2
## `geom_smooth()` using formula = 'y ~ x'

#Shift in Hatching (Days) Per 1 deg rise in temperature
  ggeffects::ggpredict(glm(hatch~Temp,data=hatch_resid_plotdata),terms=list(Temp=1))
## # Predicted values of hatch
## 
## Temp | Predicted |       95% CI
## -------------------------------
##    1 |     -6.24 | -9.10, -3.38

Graphs

Tests

hatchdate_long<- hatchdate_df %>% dplyr::select(hatch,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")

  hatchdate_corr<-hatchdate_long %>% group_by(variable) %>% cor_test(traitval,hatch,method="spearman")
  hatchdate_corr<-mutate(hatchdate_corr,p.adj=p.adjust(p,"BH"))
  hatchdate_corr$Variable<-label_df$new[match(hatchdate_corr$variable,label_df$old)]
  hatchdate_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% arrange(p.adj)
## # A tibble: 10 × 5
##    Variable              cor statistic        p    p.adj
##    <chr>               <dbl>     <dbl>    <dbl>    <dbl>
##  1 Grassland Area       0.74   106183. 9.12e-25 9.12e-24
##  2 Urban Area          -0.69   857599. 1.20e-21 6   e-21
##  3 Arable Area         -0.68   687733. 1.91e-19 6.37e-19
##  4 Minimum Temperature -0.61   784169. 6.98e-16 1.74e-15
##  5 Air Frost            0.54   224812. 3.92e-12 7.84e-12
##  6 Maximum Temperature -0.51   737671. 5.36e-11 8.93e-11
##  7 Forest Area         -0.52   624881. 6.97e-11 9.96e-11
##  8 Other Area          -0.4    713472. 4.59e- 7 5.74e- 7
##  9 Freshwater Area     -0.4    711331. 6.19e- 7 6.88e- 7
## 10 Rain                 0.18   398811. 2.99e- 2 2.99e- 2

Graphs

#Arable area
  hatchdate_plot1 <- ggplot(hatchdate_df, aes(arable_area/1000000, hatch)) +
    xlab(bquote("Arable area" (km^2))) + ylab("Hatch Date (Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.68 ***" ,x=7.5,y=125,size=7)
  hatchdate_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

#T MAX
  hatchdate_plot2 <- ggplot(hatchdate_df, aes(tmax, hatch)) +
    xlab("Mean Maximum Spawning Temperature") +
    ylab("Hatch Date (Days 1st Jan)")  +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2")+ plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.51 ***" ,x=8,y=50,size=7)
  hatchdate_plot2 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Grassland area
  hatchdate_plot3 <- ggplot(hatchdate_df, aes(grassland_area/1000000, hatch)) +
    xlab(bquote("Grassland area" (km^2))) +  ylab("Hatch Date (Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.74 ***" ,x=9,y=120,size=7)
  hatchdate_plot3
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Air Frost  
    hatchdate_plot4 <- ggplot(hatchdate_df, aes(af, hatch)) +
    xlab("Air Frost Days per Month") +
    ylab("Hatch date (Days 1st Jan)")  +
  geom_smooth(method = "lm",formula=y~poly(x,3)) + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15)  + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.54 ***" ,x=10,y=50,size=7)
  hatchdate_plot4
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

  #Urban area
  hatchdate_plot5 <- ggplot(hatchdate_df, aes(urban_area/1000000, hatch)) +
    xlab(bquote("Urban area" (km^2))) +  ylab("Hatch Date (Days 1st Jan)") +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") 
  hatchdate_plot5
## `geom_smooth()` using formula = 'y ~ x'

Individual Predictors

#Save PCA Plot
  ggsave2('Hatch PCA PLot.pdf',hatch_rotationplot,width=10,height=10,units="cm")

#Individual Predictors
  hatch_combinedplot<-plot_grid(hatchdate_plot1,hatchdate_plot3,hatchdate_plot2,hatchdate_plot4,labels=c("C","D","E","F"),nrow=2,label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Stich On Site Specific Temp Trends
  hatchdate_combinedplot2<-plot_grid(hatch_site_plot2,hatch_plot2,ncol=2,labels = c("A","B"),label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#Stich On Bayesian Model   
  hatchdate_combinedplot3<-plot_grid(hatchdate_combinedplot2,hatch_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
  ggsave2('Fig_3_Hatch Grid Plot.pdf',hatchdate_combinedplot3,width=30,height=45,units="cm")

PERCENT DEAD

Mean Trend Over Time

percdead_df %>% group_by(year,location) %>% reframe(percdead=mean(percdead,na.rm=T)) 
## # A tibble: 119 × 3
##     year location percdead
##    <int> <fct>       <dbl>
##  1  1994 T01        32.9  
##  2  1994 T02        15    
##  3  1994 T03         0.167
##  4  1994 T04         4.40 
##  5  1994 T05         2.14 
##  6  1994 T06         0    
##  7  1995 T01        33.9  
##  8  1995 T02         2    
##  9  1995 T03         2.08 
## 10  1995 T04        17.0  
## # ℹ 109 more rows
percdead_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) 
## # A tibble: 7 × 3
##   location   min   max
##   <fct>    <int> <int>
## 1 T01       1994  2010
## 2 T02       1994  2014
## 3 T03       1994  2006
## 4 T04       1994  2015
## 5 T05       1994  2014
## 6 T06       1994  2010
## 7 T11       1995  2011
hatchdate_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) 
## # A tibble: 7 × 3
##   location   min   max
##   <fct>    <int> <int>
## 1 T01       1995  2005
## 2 T02       1994  2014
## 3 T03       1994  2006
## 4 T04       1994  2015
## 5 T05       1994  2015
## 6 T06       1994  2010
## 7 T11       1996  2012
spawndate_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) 
## # A tibble: 7 × 3
##   location   min   max
##   <fct>    <int> <int>
## 1 T01       1994  2010
## 2 T02       1994  2014
## 3 T03       1994  2011
## 4 T04       1994  2015
## 5 T05       1994  2015
## 6 T06       1994  2010
## 7 T11       1996  2012
#Subset
  percdead_df_subset<-subset(percdead_df,percdead>=0)

#SampleSize 
  percdead_df %>% filter(!is.nan(percdead)) %>% nrow()
## [1] 178
# Correlations
  percdead_df %>% filter(!is.nan(percdead)) %>% group_by(location) %>% cor_test(percdead,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2     cor statistic       p conf.low conf.high method  p.adj
##   <fct>    <chr> <chr>  <dbl>     <dbl>   <dbl>    <dbl>     <dbl> <chr>   <dbl>
## 1 T01      perc… year  -0.043    -0.161 0.875   -5.27e-1     0.463 Pears… 0.882 
## 2 T02      perc… year   0.42      2.08  0.0505   2.46e-4     0.716 Pears… 0.177 
## 3 T03      perc… year  -0.34     -1.23  0.241   -7.35e-1     0.238 Pears… 0.506 
## 4 T04      perc… year   0.018     0.149 0.882   -2.23e-1     0.257 Pears… 0.882 
## 5 T05      perc… year   0.26      1.10  0.289   -2.31e-1     0.651 Pears… 0.506 
## 6 T06      perc… year  -0.17     -0.556 0.59    -6.56e-1     0.424 Pears… 0.826 
## 7 T11      perc… year   0.54      3.27  0.00301  2.09e-1     0.760 Pears… 0.0211
  # percdead_df_subset %>% filter(!is.nan(percdead)) %>% group_by(pond_id) %>% cor_test(percdead,year,method="spearman") %>% mutate(p.adj=p.adjust(p,method="BH"))
  
   
#Plot By Location 
  percdeaad_siteplot1<-ggplot(percdead_df,aes(x=year,y=percdead))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~location) + guides(fill="none") + plotopts + labs(x="Year",y="Percetage Spawn Dead") + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
  percdeaad_siteplot1
## `geom_smooth()` using formula = 'y ~ x'

#Plot By Pond ID  
  # percdeaad_siteplot2<-ggplot(percdead_df_subset,aes(x=year,y=percdead))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~pond_id) + guides(fill="none") + plotopts + labs(x="Year",y="Percetage Spawn Dead") + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
  # percdeaad_siteplot2
  
#Annotations  
  ann_text_percdead <- data.frame(year = 2001,percdead = 75,lab = "cor = 0.54 \np.adj = 0.02",
                       location = "T11")
percdeaad_siteplot1_ann<-percdeaad_siteplot1 + geom_text(data = ann_text_percdead,label = ann_text_percdead$lab,size=5)
  percdeaad_siteplot1_ann  
## `geom_smooth()` using formula = 'y ~ x'

Mean Trend with Temperature

dead_tmax_plot1<-ggplot(percdead_df)+ geom_smooth(method="lm",aes(x=af,y=percdead,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=af,y=percdead)) + guides(fill="none")  + labs(x="Mean Spawning Temperature",y="Percent Dead Spawn") + plotopts + scale_fill_brewer(palette = "Set2") 
dead_tmax_plot1
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Correlations
  spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))
## # A tibble: 7 × 10
##   location var1  var2    cor statistic      p conf.low conf.high method   p.adj
##   <fct>    <chr> <chr> <dbl>     <dbl>  <dbl>    <dbl>     <dbl> <chr>    <dbl>
## 1 T01      spawn tmax  -0.31    -1.28  0.219    -0.691    0.196  Pearson 0.307 
## 2 T02      spawn tmax  -0.45    -2.33  0.0298   -0.729   -0.0508 Pearson 0.0695
## 3 T03      spawn tmax  -0.19    -0.804 0.433    -0.594    0.288  Pearson 0.449 
## 4 T04      spawn tmax  -0.34    -2.56  0.0135   -0.561   -0.0746 Pearson 0.0695
## 5 T05      spawn tmax  -0.52    -2.52  0.0218   -0.789   -0.0892 Pearson 0.0695
## 6 T06      spawn tmax  -0.37    -1.31  0.217    -0.764    0.230  Pearson 0.307 
## 7 T11      spawn tmax  -0.15    -0.768 0.449    -0.487    0.233  Pearson 0.449

Multivariate model

#POND ID
  #with(percdead_df,table(location,LCODE))
  percdead_df$pond_id<-paste(percdead_df$location,percdead_df$LCODE,sep="_")
  nrow(percdead_df_subset) #
## [1] 178
  pdf_subtable<-table(percdead_df$pond_id)
  percdead_df_subset<-subset(percdead_df,pond_id %in% names(pdf_subtable)[pdf_subtable>9])
  nrow(percdead_df_subset)
## [1] 165
#Convert To Logits
   library(arm)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:rstatix':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loading required package: lme4
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:brms':
## 
##     ngrps
## 
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /Users/xavharrison/Library/CloudStorage/Dropbox/Mac/Documents/1. Projects/14. Frog Phenology Kat Oliver/Frog Phenology Statistics
## 
## Attaching package: 'arm'
## The following object is masked from 'package:corrplot':
## 
##     corrplot
    percdead_df_subset$percdead[which(percdead_df_subset$percdead==100)]<-99.99999
    percdead_df_subset$percdead[which(percdead_df_subset$percdead==0)]<-0.000001
    percdead_df_subset$percdead_logit<-logit(percdead_df_subset$percdead/100)
 
  
## ALL DATA, TMAX
        bf_dead<-bf(percdead_logit ~ year + ar(p=1)  +  (1|p|pond_id)) + gaussian()
        #bf_dead_perc<-bf(percdead ~ year + ar(p=1)  +  (1|p|pond_id)) + gaussian()

        #bf_af<-bf(af ~ 1 + ar(p=1) + (0+year|p|pond_id)) + gaussian()
        #bf_temp_min<-bf(tmin ~ 1 +  ar(p=1) + (1|p|pond_id)) + gaussian()
        bf_temp_max<-bf(tmax ~ year  +  ar(p=1) + (1|p|pond_id)) + gaussian()

# ### TMIN       
#   percedead_fit1 <- brm(
#   bf_dead_perc  + bf_temp_max + set_rescor(TRUE),
#   data = percdead_df_subset, chains = 4, cores = 4,
#   control = list(adapt_delta = 0.99)
# )
#     
#   summary(percedead_fit1)
#   
#   #Check 
#    pp_check(percedead_fit1,resp='percdead')
#   pp_check(percedead_fit1,resp='tmax')

### TMAX       
  percedead_fit2 <- brm(
  bf_dead  + bf_temp_max + set_rescor(TRUE),
  data = percdead_df_subset, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99)
)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
  summary(percedead_fit2)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: percdead_logit ~ year + ar(p = 1) + (1 | p | pond_id) 
##          tmax ~ year + ar(p = 1) + (1 | p | pond_id) 
##    Data: percdead_df_subset (Number of observations: 163) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Correlation Structures:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_percdeadlogit[1]     0.02      0.09    -0.15     0.19 1.00     5237     2635
## ar_tmax[1]              0.82      0.06     0.71     0.93 1.00     3348     2411
## 
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 11) 
##                                             Estimate Est.Error l-95% CI
## sd(percdeadlogit_Intercept)                     5.32      1.23     3.50
## sd(tmax_Intercept)                              1.32      0.29     0.88
## cor(percdeadlogit_Intercept,tmax_Intercept)    -0.80      0.13    -0.96
##                                             u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(percdeadlogit_Intercept)                     8.28 1.00     1457     1877
## sd(tmax_Intercept)                              2.00 1.00     1251     1869
## cor(percdeadlogit_Intercept,tmax_Intercept)    -0.45 1.00     1462     2207
## 
## Regression Coefficients:
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## percdeadlogit_Intercept   -79.58    146.65  -369.37   202.82 1.00     4398
## tmax_Intercept            -60.87     75.55  -225.05    74.43 1.00     3525
## percdeadlogit_year          0.04      0.07    -0.10     0.18 1.00     4382
## tmax_year                   0.04      0.04    -0.03     0.12 1.00     3519
##                         Tail_ESS
## percdeadlogit_Intercept     3008
## tmax_Intercept              1913
## percdeadlogit_year          2985
## tmax_year                   1913
## 
## Further Distributional Parameters:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_percdeadlogit     4.61      0.26     4.13     5.15 1.00     5159     2761
## sigma_tmax              0.54      0.03     0.48     0.61 1.00     3878     2679
## 
## Residual Correlations: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## rescor(percdeadlogit,tmax)     0.04      0.09    -0.13     0.21 1.00     5128
##                            Tail_ESS
## rescor(percdeadlogit,tmax)     3130
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Site Level Posterior Plot

percdead_fit2_site<-ranef(percedead_fit2)$pond_id

percdead_site_plotdata<- data.frame(cbind(percdead_fit2_site[,,'percdeadlogit_Intercept'][,c(1,3,4)],percdead_fit2_site[,,'tmax_Intercept'][,c(1,3,4)]))
colnames(percdead_site_plotdata)<-c("percdead","percdead2","percdead97","Temp","Temp2","Temp97")
percdead_site_plotdata$location<-substr(rownames(percdead_site_plotdata),1,3)

percdead_site_plot1<-ggplot(percdead_site_plotdata,aes(x=Temp,y=percdead))+ theme_bw() + geom_smooth(method="lm")  
percdead_site_plot2<- percdead_site_plot1 + geom_errorbar(aes(ymin=percdead2,ymax=percdead97,width=0.15),color="gray55") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=0.5),color="gray55") + geom_point(data=percdead_site_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Temperature Residual",y="Percentage Dead Spawn Residual",fill="Location") + plotopts + scale_fill_brewer(palette="Set2") + annotate(geom="text",x=-1,y=-10,label="cor = -0.8 \n [-0.96,-0.45]",size=6)
percdead_site_plot2
## `geom_smooth()` using formula = 'y ~ x'

PCA

#Removing sun
  percdead_nosun <- percdead_df %>% dplyr::select(-sun)

#Removing all rows containing NA values
  percdead_nona <- na.omit(percdead_nosun)
  nrow(percdead_nona)
## [1] 29
#Removing response variable
  percdead_noresp <- percdead_nona %>% dplyr::select(-percdead)
  
#Removing location variables
  percdead_for_pca <- percdead_noresp %>% dplyr::select(-location,-LCODE,-pond_id)
  #percdead_for_pca <- percdead_for_pca %>% dplyr::select(-LCODE)
  
#head(percdead_for_pca)
  
#Rename Variables
  percdeadpcalabs<-label_df$new[match(colnames(percdead_for_pca),label_df$old)]
  colnames(percdead_for_pca)<-percdeadpcalabs


#Performing Principal Components Analysis
  percdead_res_pca <- PCA(percdead_for_pca, graph = FALSE)
  
#Contributions
  dimdesc(percdead_res_pca, axes = 1:2)
## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Arable Area           0.9778107 7.290119e-20
## Urban Area            0.9655034 2.619852e-17
## Spawn Nitrate         0.8648117 1.450928e-09
## Maximum Temperature   0.8121883 8.833833e-08
## Minimum Temperature   0.5908527 7.387844e-04
## Rain                 -0.4438891 1.585945e-02
## Air Frost            -0.5182267 3.979574e-03
## Grassland Area       -0.9636134 5.323038e-17
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                     correlation      p.value
## Other Area            0.8275528 3.075322e-08
## Freshwater Area       0.7758503 7.625761e-07
## Minimum Temperature   0.7536917 2.358379e-06
## Year                 -0.3765206 4.409220e-02
## Forest Area          -0.3786648 4.280448e-02
## Air Frost            -0.7995917 1.958766e-07
#Producing a rotation plot
  percdead_rotationplot<-fviz_pca_var(percdead_res_pca, col.var = "black",
  
                    
  
                     repel = TRUE
  
                     ) +
    labs(x = "PC1", y = "PC2", title = "") + guides(colour="none") + plotopts

Data Processing

#Creating a data frame containing the co-ordinates of PC1 and PC2
  percdead_res_pca$ind$coord
##         Dim.1       Dim.2       Dim.3       Dim.4       Dim.5
## 1   0.7512994  1.02657507 -1.51268259  1.12409659 -1.79597090
## 2  -1.8179391 -0.62499945 -0.46850191 -1.89489783 -0.92005141
## 3  -1.6494503 -0.67708626 -0.20695827 -1.93082601 -1.19829574
## 4   0.6958885  0.18615190 -1.04519424  0.31550973 -0.78384539
## 5  -0.8130605  0.71999639 -1.78429553 -1.01966290 -0.26910012
## 6  -1.1505938  0.45776866 -1.14647909 -0.67115920 -1.12834228
## 7   1.4341075  1.10269584 -2.22329693  0.90449787  0.45993019
## 8  -1.3637829  3.92787997  0.99105864 -0.01834216 -0.01495594
## 9   3.0312263 -0.59885107 -1.29333748 -0.94173832  1.73855861
## 10 -1.5782410 -0.04932380 -1.51337000 -0.06375137 -0.62033092
## 11 -1.5306969 -0.05647557 -3.56977066 -0.71008298  2.53710954
## 12 -1.5786135 -0.04935421 -1.47285467 -0.05139918 -0.68214031
## 13  4.1095036 -0.44927316  0.43559229 -0.61179465 -0.05574561
## 14 -1.8626460 -1.75986632 -0.04224800  0.83052575 -0.19762166
## 15  4.0829075 -0.49860264  1.13086800 -0.78057490 -0.40147144
## 16 -1.0550380  3.09177794  1.50415174 -0.27695781  0.62702660
## 17 -1.0579306  3.09198789  1.69282376 -0.21865238  0.33837409
## 18  3.3708552 -1.02942816  0.97701422 -0.53933873 -0.08047223
## 19  3.9976351  0.01687874 -0.05966336  0.41406738 -0.01262993
## 20 -1.4353333  2.54179092  1.52547898  0.31284058  0.45205583
## 21 -1.4331228  2.54117462  1.50998641  0.30671756  0.47715233
## 22  3.4496714 -0.88494376  1.15779221 -0.13901098 -0.05575183
## 23  3.0210653 -1.22053408  1.13354547 -0.06159835  0.21793728
## 24 -2.0215268 -2.34070991  0.72577747  1.49505406  0.02999764
## 25 -2.9056193 -3.50928310  1.48500017  0.95991272  0.58068448
## 26 -2.2169766 -2.02808633  1.17699888 -1.13560846  0.51422729
## 27 -2.2215458 -2.02742762  1.38270306 -1.07108150  0.19851673
## 28 -1.8499598 -1.56232549  0.23147536  2.44238246  0.12439102
## 29  1.5979173  0.66189299 -0.72161393  3.03087302 -0.07923592
  percdead_pca_data <- percdead_res_pca$ind$coord[,c(1,2)]

#Creating a data frame containing the highest percentage of spawn dead per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
  percdead_glm_data <- cbind(percdead_pca_data,percdead_nona)
  #colnames(percdead_glm_data)[1:2] <- c("Dim.1", "Dim.2")
  percdead_glm_data <- as.data.frame(percdead_glm_data)
  count(percdead_glm_data) #33 data points
##    n
## 1 29

Correlations

###Model   
  percdead_glm_data %>% cor_test(.,percdead,Dim.1,method="spearman")
## # A tibble: 1 × 6
##   var1     var2    cor statistic       p method  
##   <chr>    <chr> <dbl>     <dbl>   <dbl> <chr>   
## 1 percdead Dim.1 -0.49     6065. 0.00648 Spearman
  percdead_glm_data %>% cor_test(.,percdead,Dim.2,method="spearman")
## # A tibble: 1 × 6
##   var1     var2     cor statistic     p method  
##   <chr>    <chr>  <dbl>     <dbl> <dbl> <chr>   
## 1 percdead Dim.2 -0.011     4104. 0.955 Spearman

Graphs

Tests

percdead_long<- percdead_glm_data %>% dplyr::select(percdead,spawn_nh4n,spawn_no3n,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(spawn_nh4n,spawn_no3n,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")

  percdead_corr<-percdead_long %>% group_by(variable) %>% cor_test(traitval,percdead,method="spearman")
  percdead_corr<-mutate(percdead_corr,p.adj=p.adjust(p,"BH"))
  percdead_corr$Variable<-label_df$new[match(percdead_corr$variable,label_df$old)]
  percdead_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% dplyr::arrange(p.adj)
## # A tibble: 12 × 5
##    Variable                cor statistic       p  p.adj
##    <chr>                 <dbl>     <dbl>   <dbl>  <dbl>
##  1 Arable Area         -0.53       6193. 0.00343 0.0209
##  2 Maximum Temperature -0.52       6189. 0.00349 0.0209
##  3 Grassland Area       0.49       2077. 0.00719 0.0288
##  4 Urban Area          -0.37       5571. 0.0469  0.141 
##  5 Spawn Nitrate       -0.35       5481. 0.0627  0.150 
##  6 Minimum Temperature -0.31       5323. 0.101   0.202 
##  7 Forest Area         -0.28       5192. 0.143   0.245 
##  8 Air Frost            0.26       3020. 0.18    0.27  
##  9 Other Area           0.18       3346. 0.361   0.481 
## 10 Freshwater Area      0.16       3419. 0.414   0.497 
## 11 Spawn Ammonium      -0.069      4342. 0.72    0.785 
## 12 Rain                 0.0099     4020. 0.959   0.959

Graphs

#Arable area
  percentdead_plot2 <- ggplot(percdead_df, aes(arable_area/1000000, percdead)) +
    xlab(bquote("Arable area " (km^2))) +
    ylab("Percent Dead Spawn")  +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top") + annotate(geom="text",x=7,y=75,label="cor = -0.53 \np.adj=0.02",size=6)
  percentdead_plot2
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

#TMax
  percentdead_plot3 <- ggplot(percdead_df, aes(tmax, percdead)) +
    xlab("Mean Maximum Temperature") +
    ylab("Percent Spawn \n Dead (Logits)")  +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top") + annotate(geom="text",x=7,y=75,label="cor = -0.52 \np.adj=0.02",size=6)
  percentdead_plot3 #+ facet_wrap(.~location)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

#Grassland Area  
    percentdead_plot4 <- ggplot(percdead_df, aes(grassland_area/1000000, percdead)) +
    xlab(bquote("Grassland area " (km^2))) +
    ylab("Percent Dead Spawn")  +
  geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
    theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top")+ annotate(geom="text",x=10,y=75,label="cor = 0.49 \np.adj=0.03",size=6)
  percentdead_plot4
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

####### SAVE GRID  
# percdead_grid1<-plot_grid(percdead_rotationplot,percentdead_plot2,ncol=2,labels="AUTO",label_size=25)
# percdead_grid2<-plot_grid(percdead_grid1,percdead_site_plot2,nrow=2,labels=c("","C"),label_size = 25)
# 
# ggsave2('Fig_5_Percent Dead Combined Plot.pdf',percdead_grid2,width=25,height=25,units="cm")

Individual Predictors

#Save PCA Plot
  ggsave2('Percent Dead PCA PLot.pdf',percdead_rotationplot,width=10,height=10,units="cm")

#Individual Predictors
  percdead_combinedplot<-plot_grid(percentdead_plot2,percentdead_plot3,percentdead_plot4,labels=c("C","D","E"),nrow=2,label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Stich On Site Specific Temp Trends
  percdead_combinedplot2<-plot_grid(percdeaad_siteplot1_ann,percdead_site_plot2,ncol=2,labels = c("A","B"),label_size = 25)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#Stich On Bayesian Model   
  percdead_combinedplot3<-plot_grid(percdead_combinedplot2,percdead_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
  ggsave2('Fig_4_PercentDead Grid Plot.pdf',percdead_combinedplot3,width=40,height=45,units="cm")